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pi I Abstract. We address the issue of setting up generic non-Gaussian initial conditions 

^ for N-body simulations. We consider inflationary-motivated primordial non- 

Gaussianity where the perturbations in the Bardeen potential are given by a dominant 
c/3 Gaussian part plus a non-Gaussian part specified by its bispectrum. The approach we 

I ^ I explore here is suitable for any bispectrum, i.e. it does not have to be of the so-called 

separable or factorizable form. The procedure of generating a non-Gaussian field with 
" ^ a given bispectrum (and a given power spectrum for the Gaussian component) is not 

^.f^ univocal, and care must be taken so that higher-order corrections do not leave a too 

^\ large signature on the power spectrum. This is so far a limiting factor of our approach. 

We then run N-body simulations for the most popular inflationary-motivated non- 
Gaussian shapes. The halo mass function and the non-linear power spectrum agree 
with theoretical analytical approximations proposed in the literature, even if they were 
so far developed and tested only for a particular shape (the local one) . We plan to make 
I the simulations outputs available to the community via the non-Gaussian simulations 



X 



comparison project web site http : //ice . ub . edu/-liciaverde/NGSCP . html[ 



1. Introduction 



The leading theory for the origin of primordial perturbations is inflation. Along with 
the shape and amplitude of the primordial power spectrum and the signature of a 
stochastic background of gravity waves, primordial non-Gaussianity offers a window to 
probe inflation. Even the simplest inflationary models predict deviations from Gaussian 
initial conditions, arising from the (self)-interaction of the field during inflation; for the 
simplest, slow roll, single field model these deviations are expected to be unmeasurably 
small [H E]. For a thorough review of inflationary non-Gaussianity see e.g., [3] and 
references therein. 
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An important development of the past few years is the reahzation that a large, 
potentially detectable amount of non-Gaussianity can be produced by inflation if any of 
the conditions giving the standard, single- field, slow-roll is violated. These are: a) only 
one field is responsible for generating the primordial perturbations b ) canonical kinetic 
energy of the field c) slow roll d) the quantum field was in the adiabatic (Bunch-Davies) 
vacuum. In particular the violation of each of these conditions leaves its "signature" on 
the statistical properties of the initial perturbations (e.g.,|ll El El El El El [13 E] and 
references therein). 

Deviations from Gaussianity can be characterized by the dependence of the 
bispectrum signal| on the configuration (or shape) of the three fc- vectors in its argument. 
For example local type of non-Gaussianity yields a bispectrum that is dominated by the 
squeezed configurations (where fci -C — ^3) and is generated by models of inflation 
involving multiple fields§. This shape is in general associated with models of inflation 
where non-Gaussianity is created by non-linearities developed outside the horizon [6]. 
When non-Gaussianity is created at horizon crossing during inflation, the primordial 
bispectrum is dominated by equilateral configurations {ki ~ /c2 — ks). Non-canonical 
kinetic terms will also yield non-Gaussianities of this shape. 

The non-Gaussianities produced by the most general single-field inflation models 
with approximate shift symmetry have shapes that vary from being peaked on equilateral 
configurations to being peaked on enfolded (fci ~ — ^3/2) configurations [T3]. These 
types of non-Gaussianities, have been shown [13] to be generically well described by a 
linear combination of two bispectrum shapes, one where the bispectrum signal peaks on 
equilateral configurations and another one, orthogonal to it, called orthogonal shape. 

Finally, modifications of the initial-state of the inflaton field, leave their signatures 
in a bispectrum dominated by flattened or enfolded configurations [3, [Tl] . 

The standard observables to constrain primordial non-Gaussianity are the cosmic 
microwave background (CMB) and large-scale structure. Both CMB and large- 
scale structure can measure, in principle, the bispectrum shape dependence and thus 
discriminate the shape of non-Gaussianity. The large-scale structure bispectrum 
however, includes a large contribution form non-linear gravitational evolution and 
biasing; compared with this contribution, the primordial signal "redshifts away" during 
the Universe's evolution [151 HEl [I7j . On the other hand large-scale structure can probe 
scales non easily accessible from the CMB (e.g., see discussion in [TH [19]) and offers 
other powerful probes. One technique is based on the abundance of rare events as they 
trace the tails of the underlying distribution [20j and has received renewed attention 
e.g.,[ISl [211 1221 [231 121 [2S1 [20] and references therein). Recently, the effect of non- 
Gaussian halo bias [271 128] has been demonstrated to be an extremely promising avenue 

I There are some cases where the trispectrum may be important (when, for example, the bispectrum 
is zero) but, in general, one expects the trispectrum contribution to be sub-dominant compared to the 
bispectrum one. 

§ This is also the non-Gaussianity of standard slow-roll inflation, but, as mentioned above, in this case 
the amplitude is unmeasurably small, but see e.g., [12j . 
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both for present [29l [30l EI] and future [321 |33] data. This approach uses the fact 
that the power spectrum of density extrema (where galaxies are formed) on large scales 
increases (decreases) for a positive (negative) /nl- The signal, for these two techniques 
is given by an integral over the bispectrum. The shape-dependence is thus indirect and 
suppressed, but signatures still remain [231 [12] • 

While techniques have been explored and developed to create non-Gaussian CMB 
maps with given bispectrum and power spectrum for specific and generic shapes of 
non-Gaussianities [3ll |35], to our knowledge, cosmological N-body simulations with 
non-Gaussian initial conditions specified by a bispectrum have been so far run only for 
local type of non-Gaussianity||. 

Given the importance N-body simulations have in modeling both the non-linear 
physics and observational effects that play such a crucial role in interpreting large-scale 
structure data, here we demonstrate how to set up and run N-body simulations with 
non-Gaussian initial conditions specified by a generic bispectrum. 

This paper is organized as follows. In Sec. [2] we outline how to create a three- 
dimensional non-Gaussian field with a given bispectrum. In Sec. |3] we explicitly work 
out the expressions to use in the four most popular non-Gaussian shapes discussed above 
and detail the steps for implementation. We test the non-Gaussian initial conditions in 
Sec. m In Sec. [5] we describe our simulations and use the local benchmark for 

our implementation. In Sec. [6] we present our results. We conclude in Sec. [7j 

2. Creating a 3D non-Gaussian field with a given bispectrum 

The argument of [31] (see also [35] ) which applies to a two-dimensional field expanded in 
spherical harmonics, can be generalized to apply to a three-dimensional field transformed 
to Fourier space. Pioneering work can be found in [3H1 139], here we follow a different 
route. Let's start from 

($fc,$fc,$,3) = (27r)35^(ki + k2 + k3)S(fci, h, h) . (1) 

Here B is intended to be the bispectrum of the $ field. In this convention $ is the 
Bardeen potential (i.e — Ix the standard gravitational potential, on sub-horizon scales). 
This equation refers to some time deep in the matter-dominated era i.e., we use the /nl 
CMB convention. 

We can decompose the $ field in a Gaussian and non-Gaussian parts. In Fourier 
space this reads 

$k = $^ + $r , (2) 

and therefore 

(*g$£$r) = li^^fBih, h)5^ik, + k2 + ka) . (3) 

II Simulation were run e.g. for initial conditions [36 1 137 j . but here we are concerned with inflation- 
motivated non-Gaussianities classified by their bispectrum shape. 
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If we define 



P{k2)P{kz 



with the complex conjugate of $k, and use it in Eq. (g, we recover the identity of 
Eq. 0. 

It is important to bear in mind (as already pointed out in [33]) that this procedure 
is not unique: in other words there may be more than one — non-equivalent — expression 
for yielding the same bispectrum. It is instructive to use the local non-Gaussian 
case as a worked example. The local case is usually defined in real space as [lOl [T5| HI]: 

$(x) = <|.^(x) + /nl($^(x)2 - (<l>^(x)2)) (5) 

where denotes a Gaussian random field. In Fourier space this becomes: 

$k = + /nltt.^ [ (6) 



'(27r)3 ^ 

where the constant has been absorbed in the k = mode which is ignored. 
The bispectrum for the local non-Gaussianity is: 

Bik,, k2, ks) = 2fj,T'[Pih)P{k2) + Pik2)P{ks) + P{ki)Pik,)] . (7) 

Note that Eq. ^ is not equivalent to Eq. (|4])&(|2]). Both procedures yield $ fields with 
the same bispectrum, but the $ fields obtained are different. In fact it can be shown 
that Eq. Q&Q become equivalent to Eq. (g only if B{ki, fca, k^) — > Qhi.P{,k2)P{k^) 
rather than using Eq. g. This extra "de gree of freedom" was used in Ref [3S] to 
produce non-Gaussian fields with a given bispectrum and with a minimal non-Gaussian 
contribution to the power spectrum. We return to this point below. 



3. Special cases and initial conditions implementation 

Let's start from Eq. ([T|. The bispectrum depends on the shape of the triangle made by 
the three k vectors in its argument. The detailed dependence of the bispectrum on the 
k vectors (also called in brief bispectrum shape) can help discriminate among different 
inflationary models. For example, local shape non-Gaussianities were the first type to be 
considered [ICTl [T^ l¥T] and are a direct consequence of the non-linear relation between 
the infiaton fluctuations and the curvature perturbations that couple to matter and 
radiation. While single-field slow-roll inflation generates this shape of non-Gaussianity, 
it has been shown that the amplitude is proportional to the square of the slow-roll 
parameter, which is very small [HE]- In contrast, large local non-Gaussianities can 
be generated in e.g., curvaton models [3], where the curvature perturbation can evolve 
outside the horizon or inflationary models with multiple scalar fields. Non-canonical 
kinetic terms and higher derivative contributions to the infiaton potential can produce 
significant levels of non-Gaussianity of the equilateral type if the speed of sound in these 
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models is much smaller than the speed of light, which can be realized in certain brane 
inflation scenarios |l2l HSl HU HS] . These two shapes arise due to the interactions of the 
field driving inflation. On the other hand, modified initial state (i.e. not starting from a 
Bunch-Davies adiabatic vacuum) also introduce deviations from Gaussianity. These are 
characterized by a different shape [TJ [H] which we call enfolded. Finally \13\ introduced 
the orthogonal shape: the space of non-Gaussianities produced by the most general 
single-field models, where the infiaton fluctuations have an approximate shift symmetry, 
is spanned by linear combinations of two independent shapes: the equilateral and the 
orthogonal. We shall show below that these four bispectrum types can be obtained 
from a linear combination of just three "kernels", which we will call F^"^^^, F^ and F^, 
speeding up our calculations. 

The expressions for the bispectra for the four most popular shapes mentioned above 
are reported below. 

In the local case: 

Bih, h) = 2/;^°-'F'°-'(A;i, h) ^ (8) 

where 

F'°'^-\h, k2, h) = P{h)P{h) + P{k2)P{h) + Pih)P{h) , (9) 

and we have used B to denote the bispectrum for /nl = 1, and the dependence on the 
three k vectors is understood. 
In the equilateral case: 

Bik,, k2, ks) = Qf;:iF^\k,, k2, ks) ^ f^lB^'^ (10) 

where 

= - P{k,)P{k2) + 2cyc. - 2[P{k,)P{k2)P{k,)f'^ (11) 
+ P^'''{ki)P^/\k2)P{k^) + 5cyc. 

= _ ^local _2P^ pB 

In the last equality we have introduced 

F^(fci, A;2, k^) = [P{kr)P{k2)P{k,)f' (12) 
F''{k,,k2,k,) = {[Pik^)f'[Pik2)f'Pik,) + 5cyc.} (13) 

and, for simplifying the notation, we have avoided writing down the explicit dependence 
of F on the three k vectors. 

For the template for factorized enfolded |14j : 

B{ku k2, k,) = 6f^fF^-\k,, k2, k,) ^ /^f S'^"^ (14) 



where 



^enfl _ ^local ^ ^[p(k{)P{k2)P{k:i)f/^ (15) 

- {[P{k,)Y'^[P{k2)f''P{k,) + ^cyc.] 

jplocal _|_ ^jr^ 
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Note that this is not the true enfolded configuration, but it is a factorizable template 
that approximates it well. 

Finally, for the orthogonal shape: 

B{h, h, h) = k2, h) ^ /^f (16) 

where 

F°''^{ki, k2, ks) = -3F^°" - + 3F^ . (17) 

It is important to note that these four shapes are not independent. For example 
the orthogonal one can be obtained from a combination of the other two: 

^orth p^l 2^penfi , ^orth ^^^^ 2,B^^^ (IS) 

The factorization introduced here speeds up the implementation if the initial 
conditions for N-body simulations of these four shapes: initial conditions can be 
produced directly in Fourier space and split in a Gaussian realization and three 
different non-Gaussian realizations one for the local case, one corresponding to 

F^ and one to F^ . From this one can build initial conditions for different shapes and 
different /nl- 

In the following we give some details on how we generate the initial particle 
distribution. We start off from the publicly available code by Sirko [H] for Gaussian 
initial conditions and modify it appropriately. The difference between setting up 
Gaussian and non-Gaussian initial conditions is the extra non-Gaussian contribution, 
to the gravitational potential. For a generic, non-factorizable bispectrum, we can 
compute this field by disctretizing Eq. Q, i.e., for each grid point k in Fourier space, 
we loop over all grid points k': 



r^^iJ(M-Mk.k'|)^^lg!^. (19, 



6- 

k' 

where is a random realization of a Gaussian field with the power spectrum given by 
P{k) oc fc"""^ and n<j denotes the spectral index. If the number of grid points is given by 
Ng, the computational costs of the generation of such a generic non-Gaussian field 
scales as A'^. For a modest grid size of 256^, this results already in 3 x 10^"^ evaluations 
of the summand in Eq. [19} which take about 1000 hours on a single core of a present-day 
CPU. For example, the computation for a 512^ grid would take approximately 10 days 
on 256 cores. 

For factorizable shapes the process can be greatly sped up. In fact note that in this 
case Eq. ^ can be written as a convolution (or a sum of convolutions) of two auxiliary 
fields which evaluation can be swiftly done resorting to Fourier transforms. In fact if the 
bispectrum can be written in a factorizable form as B{ki, ^2, k^) = J2i ^1(^1)^2(^2)^3(^3) 
then Eq. ^ becomes: 

1 r d^k 

^1(^)7 (^^^(k2)Q^(k + k2) (20) 



NG 
k 
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where G'{k) = bi{k)^f/P{k) and Q*(k) = b3{k)^^_^^J P{\k + k2\). The integral can 
then be quickly performed as a multiplication in real space (note the similarity with 
Eqs. (|5| and ^ where instead of the $ field we have two auxiliary fields G and Q). 



While for factorizable shapes the direct summation approach (Eq. 19) and the 
convolution approach are mathematically equivalent in a practical implementation they 
may be affected by different numerical effects. Here we explore both implementations 
for the local type. In the other cases we implement the direct summation approach to 
explore whether this approach is viable and uncover possible bottlenecks or limitations. 

Next, the linear density field 5k at 2; = is derived from the potential •I'k through 
the transfer function and the Poisson equation: 

^ _ 2k^T{k)D{z)^ ^21) 
3 VL^Hq 

where D{z) is the linear growth function normalized to (1 + ^) in the matter-dominated 
era, and T{k) is the transfer function obtained with CAMB [17] and normalized to unity 
on large scales. VLm is the present-day matter fraction and Hq the Hubble constant. The 
particles are then displaced from a regular grid according to the displacement field at 
the initial redshift, Zi = 49, using the Zel'dovich Approximation^. 

4. Testing the initial conditions 

In this section we analyze the quality of the non-Gaussian initial conditions. First 
we consider the 1-point function, i.e. the probability distribution function of the density 
contrast 5 at the position x. Especially, we compute the variance, (6^), and the skewness, 
(5^), as a function of the smoothing scale R and compare the results with the analytic 
predictions. Furthermore, we calculate the power spectrum of the non-Gaussian initial 
conditions and demonstrate that the deviations from the Gaussian case are in most 
cases small. 

First we generate eight Gaussian realizations of the density field on a grid of size 
256^ in a box with a side of 600 Mpc/h. From these Gaussian realizations we produce for 
each of the previous mentioned types of non-Gaussianity, i.e. local, equilateral, enfolded, 
and orthogonal, a non-Gaussian density field with an /nl of —500, —250, —100, 100, 
250, and 500. 

For the local type, we compute the non-Gaussian contribution in three different 
ways: a) the traditional way by squaring the Gaussian density field in real space b) in 
Fourier space using our ansatz Eq. (19) with the bispectrum given by Eq. ([T]) using 



again Eq. (19) but with B{ki, k2, k-^) — > 6/nl-P(^2)-P(^3)5 this recovers the traditional 
method in real space except for aliasing effects introduced by the finite grid size. 

In Fig. |4] and Fig. [2] we show the variance and skewness of each realization 
(data points) for all types of non-Gaussianity with /nl = 100 and compare them to 

% Since we will be interested in the non-Gaussian to Gaussian ratio of our statistics, the detailed 
implementation of the initial displacement field (i.e., Zel'dovich or Second-Order Lagrangian 
perturbation theory) does not matter. 
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Figure 1. Variance of the linear density field at z — as a function of the smoothing 
radius R. The non-linearity parameter is /nl = 100. The symbols show the variance 
derived from the 8 different realizations of each type of non-Gaussianity. The lines 
depict the theoretical predictions. For clarity, the case "local a)" and "local b)" are 
slightly shifted horizontally. All other symbols and lines fall on top of each other. 



analytic predictions (solid lines). The analytic prediction for the skewness is obtained 
by integrating the bispectrum (see e.g., Eq. (4.13) of Ref. 118] ). We have truncated 
the integrals at the maximum and minimum k sampled by the simulation box. The 
magnitude of the effect depends on box size, scale and type of non-Gaussianity. For our 
simulation settings, we find it always to be below 15%. 

The moments of the density are computed from the linear density field smoothed 
with a top-hat filter of radius R. In case of the skewness, the (6^) of the Gaussian 
realization is subtracted from the total skewness in order to reduce the noise introduced 
by the finite volume and grid size. The variances agree well with the theory except for the 
type "local b)". The increased variance for this type is caused by the term ($^*^$^^), 
which, in this case, is not negligible. This term gives rise to a P^(fc) term multiplied 
with a divergent integral, which in our case of a finite box and grid is truncated. This is 
reminiscent of the discussion in sec. V of [35]; for more details see Appendix A Apart 
from the "local b)" case, the values of the skewness obtained from the initial conditions 
also agree well with the predictions. Only for the orthogonal case small deviations at 
larger radii are visible. 

Next we consider the power spectrum of the initial particle distributions. We 
compute the power spectrum by first assigning the particle to a 512^ grid using the 
cloud-in-cell scheme. Then we perform a Fast Fourier Transform and average the k- 
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Figure 2. Skewness of the linear density field at z = for different types of non- 
Gaussianity with /nl = 100. Symbols show the skewness of the 8 different realizations. 
Note that the skewness of the Gaussian field due to the finite volume and grid size 
is subtracted from the measured skewness of the non-Gaussian field. Lines show the 
theoretical predictions. Only the orthogonal case has a negative skewness for positive 
/nl- For clarity, the case "local a)" and "local b)" are slightly shifted horizontally. 



modes over spherical shells with a thickness of Ak = 0.01 h/Mpc. The ratios of the 
power spectrum of the non-Gaussian and the Gaussian initial conditions are shown in 
Fig. |3j The eight different realizations are depicted by the different symbols. For clarity, 
we show the individual realization only for k < 0.1/i/Mpc. The solid line represents 
the mean of the realizations. On large scales, the power spectrum of the "local b)" 
case deviates strongly from the Gaussian power spectrum. At the lowest wave number 
k = 0.01 h/Mpc the power spectrum is enhanced by almost an order of magnitude! This 
behavior explains the offset of the variance and can be traced back to the second-order 
term ($^'^$^^), which is further explored in the Appendix A We observe a similar 
effect for the orthogonal shape, although the deviations are much smaller and a possible 
systematic shift, i.e. enhancement or suppression of the power spectrum, is — if at all — 
at the few percent level. Nevertheless we discard the orthogonal shape from now on and 
refer to the Appendix A for more details on this issue. The case "local a)" and "local c)" 
agree perfectly with each other as it is expected from their mathematical equivalence. 
Both of these implementations of the local type and the equilateral case (almost not 
visible) do not show deviations from the Gaussian power spectrum. The power spectrum 
of the enfolded shape has very small deviations from the Gaussian power spectrum, the 
mean only deviates at the sub-percent level. We also checked variance, skewness, and 
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power spectrum for the other /nl values and found quahtatively similar results. The 
deviations from the Gaussian power spectrum scale roughly linear with /nl, except for 
the "local b" case for which they scale as 



see 



Appendix A for details). 




Figure 3. Ratio of the power spectra derived from the particle distributions of 
the non-Gaussian (/nl = 100) and Gaussian initial conditions. For k < 0.1/i/Mpc, 
individual realizations are shown by symbols. The lines depict the arithmetic mean 
of the realizations. For clarity, the case "local a)" and "local b)" are slightly shifted 
horizontally. 



5. Testing the simulation settings 

In this section we perform several tests to assure that the settings — like force and 
mass resolution, box size, and initial redshift — of our simulations are good enough to 
assure that systematic effects are smaller than the statistical errors. Note that we only 
consider ratios of non-Gaussian to Gaussian quantities in this paper. In particular, we 
are interested in ratios of the non-linear power spectrum and the halos mass function. 
Most of the systematic effects due to numerical limitations are expected to cancel out 
when considering ratios (see e.g. |19]). In addition, since the non-Gaussian and Gaussian 
simulations are based on the same realization of the random density field, also the 
statistical error on the ratios is reduced vastly. 

The standard settings of our main set of simulations can be found in the first half of 
Tab. [1} The second half of the table describes the simulations which we use to investigate 
the numerical uncertainties of our main set of simulations. The generation of initial 
conditions done with the method "local a)" — the traditional real-space implementation 
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of local non-Gaussian initial conditions — is by far not as computationally intensive as 
the method in Fourier space described in Sec. |2]and|3j We exploit this fact and use the 
case "local a)" to perform larger N-body simulations to explore the numerical limitations 
of the smaller simulations of the main set. In addition, we use the case "local a^ as a 
benchmark for the implementation of local non-Gaussianity in Fourier space. 

All our N-body simulations are performed with the publicly available code Gadget- 
2 [50] . Our choice of cosmology is a flat ACDM cosmology, which is consistent with the 
seven-year WMAP results [ST]. In particular, we choose the following values: = 0.27, 
flbh^ = 0.023, h = 0.7, = 0.95, and as = 0.8, where Qm is the matter density, Qb the 
baryon density, h the Hubble parameter, Us the spectral index of the primordial power 
spectrum, and as the rms of linear density fluctuations at z = in a sphere of 8 Mpc/h. 

First, we consider the ratio of the non-linear power spectra derived from the local 
non-Gaussian simulations with /nl = 100 and Gaussian simulations at z = 1 and z = 0. 
The results are presented in Fig. |4j The black dashed and sohd lines depict the average 
of the eight realizations of the case "local c)" at redshift z = 1 and z = 0, respectively. 
The magenta and cyan lines show the two different implementations of the local type of 
non-Gaussianity for a single realization. The agreement is excellent and demonstrates 
the functionality of our method in Fourier space. The ratios derived from simulations 
of a larger box (1200 Mpc/h, blue lines) and from simulations with a higher starting 
redshift {zi = 99, red lines) are consistent with the ones derived from our main set of 
simulations. Only in the case of higher force and mass resolution, we find statistically 
significant deviations on smaller scales, i.e. the ratio falls off for k > 0.7 and k > 0.3 at 
z = 1 and z = 0, respectively. The small reduction of the ratio is probably caused by 
the enhancement of non-linear power due to the better resolution. 

The second quantity we are interested in is the ratio of the non-Gaussian and 
Gaussian halo mass function. To identify halos in the simulation data, we use the 
publicly available halo finder AHF [521 ES]- AHF defines halos to be gravitationally 
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Figure 4. Ratio of the power spectrum of local non-Gaussian (/nl ~ 100) and 
Gaussian simulations at z = 1 (dashed lines) and z = (solid lines) for different 
simulation settings (see Tab.[T]). The dashed and solid black lines show the mean ratio 
derived from the main set of our simulations of type "local c)". 

bound objects which have a spherical overdensity (SO) given by the redshift dependent 
virial overdensity. The minimal number of particles in a halo is 20, thus in the main 
set of our simulations we find halos of mass 2 x 10^3 MqIK or higher. Ahhough 20 
particles are too few to reliably resolve all halos of the corresponding mass, we find 
that the ratio of the mass functions is not affected by the low mass resolution. This is 
demonstrated in Fig. [5] The black dashed and solid lines show the averaged fractional 
difference in the non-Gaussian and Gaussian mass functions at redshift 2; = 1.5 and 
z = 0.67, respectively, derived from the eight simulations of the type "local c)". The 
results obtained from the simulations with an eight times higher mass resolution are 
depicted by the green symbols and are consistent with the lower resolution results. In 
addition. Fig. [sjshows that the ratio is not biased by the finite-volume (blue symbols) nor 
the initial redshift (red symbols). The very good agreement between the two different 
ways of setting up the local case (cyan and magenta symbols) reassures us of the correct 
implementation of the method described in Sec. [2] and [3| 

Overall the results of this section give us confidence that for our main set of 
simulations — including the other types of non-Gaussianity — systematic effects due 
to numerical limitations are within the statistical errors. 



z=1, higher resolution 





Figure 5. Fractional difference in the mass function of local non-Gaussian (/nl = 100) 
and Gaussian simulations at z = 0.67 (filled circles) and z = 1.5 (open squares) for 
different simulation settings (see text). Error bars represent Poisson errors. The dashed 
and solid black lines show the average of the fractional difference derived from the main 
set of our simulations of type "local c)". 



6. Results 

Here, we present our findings derived from the non-Gaussian simulations of the local, 
equilateral, and enfolded type for /nl = -500, -250, -100, 100, 250, and 500. First 
we present the results for the non-linear matter power spectrum. Afterwards, we turn 
to the halos mass function. The investigation of the non-Gaussian halo bias effect e.g., 
[27| [28] and the bispectrum is left to forthcoming work [16] . 

6.1. Non-linear power spectrum 

High precision in the theoretical prediction of the nonlinear matter power spectrum at 
~ 1 h/Mpc and above is needed to derive unbiased results from the data of upcoming 
weak-lensing surveys e.g., [5l]- Future Lyman-alpha forest surveys will also access these 
scales and test them with small statistical error-bars. Non-Gaussianities in the initial 
conditions alter the nonlinear power spectrum at the few percent level. In Fig. |6] (local). 
Fig. [T] (equilateral), and Fig. [s] (enfolded), we show the ratio of the non-Gaussian and 
Gaussian power spectrum for different type of non-Gaussianities with /nl = 500 (green), 
250 (red), and 100 (blue) at redshift z = 1 (open squares) and z = (filled circles). 
The black dashed lines correspond to the perturbation theory prediction using the time- 
renormalization group (TRG) approach |55]. For z = 1 and k < 0.2 h/Mpc, the TRG 
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Figure 6. Ratio of tiie non-Gaussian and Gaussian power spectrum for the local type 
for different values of the non-linearity parameter /nl at z = 1 (open squares) and 
z — 1 (filled circles). The dotted black lines show the Time Renormalization Group 
(TRG) perturbation theory prediction of [55 . Note we scaled the /nl = 80 data of 
[55] linearly to /nl = 100. The error bars show the rms from the 8 realizations. 
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Figure 7. Same as in Fig. [6] but for the equilateral type of non-Gaussianity. 
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Figure 8. Same as in Fig. [6] but for the enfolded type of non-Gaussianity. 

predictions agree very well with the results of the N-body simulations. At z = 0, 
perturbation theory slightly overpredicts the effect on scales around k ~ O.l/i/Mpc 
before it breaks down on smaller scales {k < 0.2 h/Mpc). 

Note that the maximum of the ratio is larger at higher redshifts and its positions 
is shifted to smaller scales. The shape, especially the peak location and peak height, is 
consistent with predictions of the halo model presented in |56j . 



6.2. Halo mass function 

At the high-mass end, the halo mass function is very sensitive to a primordial skewness 
in the probability distribution function of the density field ^57j|. Hence, galaxy cluster 
surveys offer in principle the possibility to probe primordial non-Gaussianities. In Fig. |9] 
we present the ratios of the cumulative mass functions derived from the non-Gaussian 
and Gaussian simulations. We show the results for the local, equilateral, and enfolded 
shape with different /nl at redshift z = 1.5 and z = 0.67. In addition to the data points 
we plot the analytic predictions of [20] (MVJ, blue lines) and [IS] (LV, red lines). In 
order to convert the analytic ratio of the non-Gaussian and Gaussian mass functions, 
Tj^ciM), provided in [20] and [18] into the corresponding ratio of the cumulative mass 
functions, Rng{> M), we use the fitting formula of [58] for the SO halo mass function, 
nTinker(M), i.e. 



R 



NG 



>M) 



/m ^Tmker(M)rfM 



(22) 
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Figure 9. Ratio of tlie non-Gaussian and Gaussian cumulative mass functions for 
the local, equilateral and enfolded shape at redshift z — 1.5 and z — 0.67 for different 
values of /nl- The symbols depict the results derived from the simulations using 8 
different realizations for each type of non-Gaussianity. For the error bars, we assume 
Poisson errors. The blue solid lines show the predictions of [20l and the red lines 
represent the model of [18]. 



We checked that the integrated Tinker fit is a good fit to the cumulative mass function 
of the SO halos found in our Gaussian simulations. 

Note that, compared to [21], we do not find that substituting the linear spherical 
collapse overdensity by 6c — > y/<j6c (with q = 0.75) improves the agreement between the 
N-body data and the predictions. This is in accordance with the findings of [SHI EO] , who 
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argue that the differences are due to the different halo identification algorithms. While 
in this paper we use halos defined by spherical overdensities, [21] applied a Friends-of- 
Friends halo finder to their simulations. 

At the high-mass end, we find that LV fits the data better for positive /nl, whereas 
MVJ gives better predictions for negative /nl- In general, for large positive /nl the 
theoretical predictions overestimate the increase of halo number density at the low- 
mass end. However, for all shapes, the qualitative behavior is captured well by the 
analytic models. There is an indication that using spherical overdensities halos a value 
q ~ 0.9 yield a slightly better fit to the simulations data. Given our simulations size, 
this effect is however not highly significant and only visible for high values of /nl- While 
the exact value of q is an important theoretical issue, recall that a 10% uncertainty in q 
will propagate into a 10% error on /nl- For values of /nl ^ 10, as expected from most 
models, this uncertainty is at most comparable to the expected statistical errors. 

7. Conclusions and discussion 

We have addressed the issue of setting up generic non-Gaussian initial conditions for N- 
body simulations. In the era of precision cosmology, N-body simulations have become a 
crucial tool to test, develop and calibrate any statistical analysis of large-scale structure 
surveys. While the approach of constraining primordial non-Gaussianity with large-scale 
structure and galaxy surveys has recently received renewed attention, as far as we know, 
N-body simulations have been run only for the so-called local-type of non-Gaussianity. 
The shape of non-Gaussianity, however, is crucial if one wants to use such a signal as a 
window into the generation mechanism for primordial perturbations. 

Building on the expertise developed in the context of Cosmic Microwave 
Background non-Gaussianity, we have shown how to set up non-Gaussian initial 
conditions for any type of non-Gaussianity specified by a primordial bispectrum. Given 
the current cosmological constraints all non-Gaussian fields we consider are given by the 
sum of a (dominant) Gaussian component and a (subdominant) non-Gaussian one. 

The implementation is based on direct summation in Fourier space and is 
significantly more computationally intensive than the workhorse local case defined in 
real space, as it scales as where iV| is the number of grid points in the simulation 
box. We have investigated the numerical effects of such an approach by comparing 
the results of the local-type case obtained with the two approaches. For factorizable 
templates, the implementation can be sped up significantly: in fact the operation can 
be rewritten as a convolution of two suitably defined auxiliary fields. This is useful, 
however not the main goal of our investigation, as our aim was to develop an approach 
suitable for non-factorizable templates. This is relevant because for some physically 
motivated inflationary models the existing factorizable templates may not be a good 
approximation of the effect e.g., [631 EH EH] and references therein. 

The procedure of generating a non-Gaussian field with a given bispectrum (and a 
given power spectrum for the Gaussian component) is not univocal and care must be 
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taken so that higher-order corrections do not leave a too large signature on the power 
spectrum especially on large scales. This is so far a limiting factor of our approach, as 
these components can be kept under control only in specific cases. More investigation 
is clearly needed, possibly along the lines proposed by [35] . 

We have then explored the effects of several popular forms of primordial non- 
Gaussianity on the halo-mass function and the non-linear power spectrum deferring the 
analysis of other statistics such as the non-Gaussian halo bias to a forthcoming paper. 
We confirm that the non-Gaussian correction to the halo mass function is determined 
by the primordial skewness and that a suitable combination of the different analytical 
approximations proposed in the literature (depending on the regime of applicability) 
offer a good fit to the simulations data. We also confirm that independently of the type 
of non-Gaussianity spherical-overdensity halos need a q correction factor much closer to 
unity (if any correction at all) than Friends-of-Friends halos. 

We believe that the methodology illustrated and developed here will be relevant 
for the on-going and planned efforts of constraining primordial non-Gaussianity from 
large-scale structure including surveys of galaxies, high-redshift clusters, Lyman-alpha 
forest and weak lensing. 
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Appendix A. Non-Gaussian contributions to the initial power spectrum 

As before we write the gravitational potential as the sum of a Gaussian and non-Gaussian 
field, $k = + ^k^*^) where the non-Gaussian part is quadratic in The power 
spectrum, defined by ($k'^'q) = (27r)'^5-^(k + q)P{k), consists then of the following 
terms _ /nl($^$^$^), and ($^g$jvg^ _ /2L($«$«$f?$^). The 

term (<|)'^<1)^'^) will vanish in theory as it involves an odd number of Gaussian fields. 
In practice, we have only a finite number of modes to perform the average. Hence, 
especially on larges scales where there are only a few modes in the box, this term is not 
exactly zero. We investigate this term further below. 

Since the term (<|)^'^<|)^'^) is of second order and the potential is small, one could 
think that this term is negligible for reasonable values of /nl- However, our formula for 
^NG involves integrals which — depending on the bispectrum — are divergent. In order 
to explore this further, we derive the general expression for ($^'^$^'^). After applying 
Eq. Q for (^^'^^ using the definition of the power spectrum, and integrating over the 
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delta function we obtain 

As a first example let us consider the case in which the bispectrum of the local type 
(Eq. |8]) is used for B in Eq. (|4]). In this case, we get three terms 

= ^/NL^^(k + q) {/ d'k'P{k')P{\k + k'l) 
+ 4P(A;) J d^k'P{k') 

+2P^(fc)|d3^'(^il^ + l)}, (A.2) 



where the first and second term are very similar to the ones discussed in 102] ; the 
second term is just a fc-independent renormalization, which for reasonable values of /nl 
is negligible small because of the truncation of the integral due to the finite volume of the 
simulation box, while the first term gives rise to a fc-dependent renormalization which 
results in a change in the slope of the power spectrum. However, for realistic values of /nl 
this change in the slope is very small [62]. The last term, proportional to P'^{k), causes 
the largest modification of the power spectrum. Even for still allowed values of /nl this 



term gives significant contributions to the power spectrum on large scales (see Fig. Al ). 
In order to circumvent this problem, we consider B{ki, k2, k^) — y Qf^-LP{k2)P{kz)-, 
which to first order produces the same bispectrum, but does not imply this kind of 
divergent integral. For this ansatz we obtain only the aforementioned scale-dependent 
term 

i^^NG^NG^ = 2fl^5''{k + q) 1 d'k'P{k')P{\k + k'l) . (A.3) 



In Fig. Al , we show the power spectrum for an extreme /nl to make the change in the 
slope visible. 

Now we turn to a more general case consisting of combinations of pi°=^' (Eq. [9|, 



and F (Eq. 12). To do this we calculate the different products of these three functions. 



The arguments of the F-f unctions are k, k\ and k = |k + k'| 



P{k')P{K) 
P{k')P{K) 



P^'\k) [P^'\k')P^''\K)\ (A.4) 

P{k) [P^'\k')P^'\K) + P^/\k')P^'\K)\ 
+ P^/\k) [p2/3(^)+p2/3(^/)] 

+ P'''\k) [P^/\k) + P^'\k')] (A.5) 
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Figure Al. The power spectrum of the primordial potential obtained from a single 



realization of $k for different choices of B used in the ansatz Eq. (19 1. The potential 
field $!(• was realized on a grid of size 256"^ in a (600Mpc//i)^ box. The dotted black 
lines show the theoretical predictions. For "local b)" and "local c)", the deviations 
from the Gaussian case originate from the term ($^'^$^'5), while for the orthogonal 
case the deviations are mainly caused by the term ($'-^$^''). 



A IT'local 



P{k')P{K) 



+ P'''\k) [p-^l\k')P''/\K) + p2/3(^')p-l/3(^)] (A.6) 



P{k')P{K) 



P'"\k) [P^'\k')P{K) + P{k')P^I\K) + 2p2/3(A;')p2/3( 

+ P{k) [2P{k') + 2P{n) + 2p2/3(fc')pi/3(^) + 2P^/\k')P^'\K 

+ P^'\k) [p-^l\k')P{K) + P{k')p-^''\K) 

+2P^'\k') + 2p2/3(^) + Qpy\k')P^'\K 
+ P''l\k) [2P-^'\k')P''l\K) + 2P^'\k')p-^l\n) 

+2P^'\k) + 2P^/\k')\ 
+ P\k) [P^I\k')P~^'\n) + p-^l\k')P^'\K) + 2] 



(A.7) 
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+ P{k) [P^/\k')P^/\K) + P^'\k')P^/\K)\ 

+ P^/\k) [p-^/\k')P{K) + P{k')p-^'\K) 
+p2/3(fc')+p2/3(^)j 

+ P^/\k) [p-^/\k')P{K) + P{k')p-^'\K) 

+py\K)+P^'\k')\ 
+ P\k) [P^'\k')p-^'\K) + p-^'\k')P^'\K) 

+P'"\k')P-^'\K) + p-^'\k')P''l'\K)\ (A.8) 



;ilocal TTilocal 



P{k')P{K) 



[Pik')PiK)] 

+ P{k) [2P{k') + 2P{K)] 

+ P\k) [P{k')p-\K) + p-\k')P{K) + 2] 

(A.9) 

In the limit of <^ fc', i.e. k k', all these second-order terms can be written 
as P'^{k)P'^~^{k') with < r < 2. Using P{k) ^ k~^ and performing the truncated 
integration over k' , we find that on large scales the contribution of the terms increase 
significantly with larger r. 

For certain combinations of the F-functions many of the above terms cancel. For 
example, with B oc (^F^°^^^ — F^) the P^(/c)-terms vanish in the above limit. Moreover 



for the equilateral case (see Eq. 11) all terms with r > 1 cancel, whereas for the enfolded 
type the term with P'^^^(k) does not vanish. 

Now let us return to the term (<|)'^<1)^'^) linear in /nl- In the discretized version 
this becomes: 

For small k, there are only a few modes in the box and the skewness of the realized 
Gaussian field is not completely vanishing. Depending on the choice for B this noise can 



be amplified significantly. In particular, for our ansatz of the orthogonal type (Eq. 17) 
this contributes substantially to the power spectrum and causes the deviations visible 
in Fig. [3j 
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